Quantum phase slips in one-dimensional superfluids in a periodic potential 
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We study the decay of superflow of a one-dimensional (ID) superfluid in the presence of a periodic 
potential. In ID, superflow at zero temperature can decay via quantum nucleation of phase slips even 
when the flow velocity is much smaller than the critical velocity predicted by mean-field theories. 
Applying the instanton method to the 0(2) quantum rotor model, we calculate the nucleation rate 
of quantum phase slips V. When the flow momentum p is small, we find that the nucleation rate 
per unit length increases algebraically with p as T/L oc p 2K ~ 2 , where L is the system size and K 
is the Tomonaga-Luttinger parameter. Based on the relation between the nucleation rate and the 
quantum superfluid-insulator transition, we present a unified explanation on the scaling formulae 
of the nucleation rate for periodic, disorder, and single- barrier potentials. Using the time-evolving 
block decimation method, we compute the exact quantum dynamics of the superflow decay in the 
ID Bose-Hubbard model at unit filling. From the numerical analyses, we show that the scaling 
formula is valid for the case of the Bose-Hubbard model, which can quantitatively describe Bose 
gases in optical lattices. 

PACS numbers: 03.65.Xp,03.75.Kk, 03.75.Lm 



I. INTRODUCTION 



Recently, superfluidity and superconductivity in one 
dimension (ID) have been experimentally studied in var- 
ious systems, including superconducting nanowires [il- 
liquid helium in nanoporcs 
atoms in optical lattices 




and ultracold bosonic 
A common property 



found in these different systems is that the transport in 
ID is significantly suppressed compared to that in higher 
dimensions. This suppression of the transport might be 
interpreted as a consequence of stronger effects of thermal 
and quantum fluctuations in ID. At temperatures higher 
than a certain characteristic value, thermal fluctuations 
allow the amplitude of the superfluid order parameter 
to vanish and its phase to unwind, leading to the de- 
cay of superflow [D3 - [l4| . Such a process is often referred 
to as a phase slip. When the temperature is sufficiently 
low, thermal fluctuations are suppressed and the nucle- 
ation of phase slips due to quantum tunneling provides 
dominant contributions to the superflow decay [l5l - |23| . 
Indeed some of the experiments with superconducting 
nanowires have observed the crossover from the regime 
of the thermal activation to the quantum regime |l|-|5|. 
Moreover, the experiments studying the transport of ID 
Bose gases in optical lattices @ showed a good agreement 
with the theoretical analyses at zero temperature by the 
time-evolving block decimation (TEBD) method [13, HE] 
that accurately takes into account the effects of quantum 
fluctuations. In light of this agreement, it is highly likely 
that the quantum regime has been achieved also in cold 
atom systems. Therefore, it is important to accurately 
calculate the nucleation rate of quantum phase slips for 
ID superfluid. 

In particular, the decay of ID superflow in a periodic 
potential have attracted much interest because the trans- 
port of ID atomic Bose gases has been studied in the 



presence of an optical lattice. In addition, in the exper- 
iments of liquid helium absorbed in nanopores, an inert 
layer of solid helium covering the wall of thepores acts as 
an external potential for ID liquid helium [7j , which may 
be regarded as a periodic potential [2(| . Having in mind 
ultracold atom experiments, most of previous theoretical 
studies have analyzed transport properties in the pres- 
ence of a parabolic trapping potential in addition to pe- 
riodic potentials by means of various numerical methods 
beyond mean-field approximations, such as exact diago- 
nalization [27j , truncated Wigner approximation [H, [2t| , 
fermionization method [30l . [3l| . and TEBD (or equiv- 
alently time-dependent density-matrix renormalization 
group) [13, HE Hl|- They have not explicitly related 
the ID transport at zero temperature to quantum phase 
slips. In contrast, in Ref. [221 ] the authors have pointed 
out the connection with phase slips and used the instan- 
ton method to analytically calculate the nucleation rate 
in a homogeneous lattice system. However, their instan- 
ton analyses in ID have been restricted to the parameter 
regions far from the Mott transition and close to the su- 
perfluid critical velocity predicted by mean-field theory. 

In the present paper, we investigate quantum phase 
slips in ID superfluids in the presence of a periodic po- 
tential. Applying the instanton techniques to the 0(2) 
quantum rotor model, we calculate the nucleation rate 
r as a function of the flow (quasi-)momentum p. Since 
we treat the phase degrees of freedom on all the sites as 
independent variables in contrast to previous work that 
uses a single-collective- variable approximation [HJ, we 
can analyze the entire region of the momentum. Espe- 
cially, when the momentum is much smaller than the 
critical value p c , we find that the nucleation rate per 
site obeys T/L oc (p/p c ) 2K ~ 2 , where L is the number 
of lattice sites and K is the Tomonaga-Luttinger pa- 
rameter. This power-law behavior with respect to the 



2 



momentum means that the lifetime of supcrflow can be 
practically infinitely long if p <C p c , namely the presence 
of superfluidity in ID. It should be noted that a similar 
power-law behavior of the nucleation rate has been pre- 
viously found in the case of ID superfluid in the presence 
of a single-barrier potential [HI, [lj| or a disorder poten- 
tial [2l[, i.e. T oc (p/p c ) 2if_1 for a single barrier and 
r/L oc (p/pc) 2K ~ 1 in the case of disorder. We emphasize 
that in the case of periodic potentials the exponent takes 
a different value. We discuss this difference from a view- 
point of the relation between the nucleation rate and the 
quantum superfluid-insulator transition. Moreover, we 
investigate the real-time dynamics of the superflow de- 
cay via quantum phase slips in the ID Bose-Hubbard 
model (BHM) at unit filling by means of the quasi-exact 
numerical method of TEBD [H|. From the TEBD 
calculations, we confirm the validity of the scaling for- 
mula in the model that is quantitatively relevant to ex- 
periments with ultracold Bose gases confined in optical 
lattices. 

The remainder of the paper is organized as follows. In 
Sec. HIl we introduce BHM describing ID bosons in a 
periodic potential and briefly review a derivation of the 
0(2) quantum rotor model from BHM. In Sec. IIII1 we 
calculate the nucleation rate of quantum phase slips us- 
ing the instanton techniques and find a scaling formula 
of the nucleation rate with respect to the momentum. In 
Sec. IIV1 we discuss a physical interpretation of the scal- 
ing formula from a viewpoint of the relation between the 
phase-slip nucleation and the superfluid-insulator transi- 
tion. In Sec. El applying TEBD to BHM at unit filling, 
we compute the quantum dynamics of superflow, which 
exponentially decays in time. We numerically extract 
the nucleation rate of quantum phase slips and compare 
it with the scaling formula. In Sec. I VII we summarize 
our results and describe possible future work. 



II. 0(2) QUANTUM ROTOR MODEL 

We consider a system of ID lattice bosons in a ring- 
shaped geometry, i.e. with a periodic boundary. To study 
such a system, one often uses BHM that especially allows 
for a quantitative description of Bose gases confined in 
optical lattices. The quantitative validity of BHM has 
been confirmed for both equilibrium and dynamical prop- 
erties through the thorough comparisons between exper- 
iments on Bose gases in optical lattices and exact numer- 
ical analyses on the Bose-Hubbard model [35U381 ] . How- 
ever, before directly analyzing BHM, in order to gain 
analytical insights we use the 0(2) quantum rotor model 
that can be derived from BHM under a certain condition. 
The formulation of the instanton method is much simpler 
for the quantum rotor model [22, [3!| . In this section, we 
review a derivation of the quantum rotor model. 

We start with the grand canonical partition function, 



Z = / Vb*Vbexp 



S[b*,b] 



(1) 



where the action S[b*, b] for the BHM is given by 
S\T,b]=£f**r 

]=1 J -— 

b*{r)h^-b {r) - J (b*(r)b j+1 (r) + &* +1 (t)^(t)) 



+ ^*{r)b*{r)b ] {T)b ] {T) - nb*{T)b 3 (T) 



(2) 



where U is the onsite interaction, J the hopping energy, 
and L the number of lattice sites. Here, fi w Uv is the 
chemical potential and v is the filling factor. For con- 
venience we introduce finite small temperature T corre- 
sponding to the inverse temperature /3 = (fcflT) -1 . In 
the end of calculations we will take the limit of T — > 0. 
Inserting bj = y /nje iej , the action is rewritten as 

S[n,9] = > ' dr hn-j 



E 

3=1 

-2y/nJnJ+iJ cos (9j+i 



dr 



1 drij 
2rij dr 



9j) + ^K" 



.(3) 



We split the number of particles per site into its aver- 
age and fluctuation as rij = 



Srij, and assume that 



Uv 3> J and v 3> Srij. Then, we find that the action is 
approximated as 



S\n, 



E 



dr 



r ^ d 3 
ih(v + bnj)—^- 

OT 



-2i/Jcos (9 



3+1 



(4) 



Since Eq. (jj) contains only the linear and quadratic terms 
with respect to number fluctuations Srij, these degrees 
of freedom can be integrated out. Then, the action is 
described in terms of the phases as 



S[6] = J2 / 2 dr 



OT 



2U 



dr 



-2vJcos{9 J+1 - 9j)} 



(5) 



When the filling factor is irrational, the first term makes 
the net contribution of the instanton or bounce solutions 
to the partition function to be zero [io[ • When the filling 
factor is integer (commensurate filling), the first term is 
necessarily equal to h x 2nl, where I is an integer, and its 
contribution to the partition function is unity regardless 
of the trajectory of 9j. In the latter case, the effective 
action takes the form of the quantum rotor model, 



S[6} = 



L 

E 



dr 



2U 



, 89, , 
I — - ) 



-Oj) 



.(6) 



Introducing the dimensionless parameter K = iry/2vJ/U 
and the sound velocity u = d\j2vJU /h, the action can 
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be rewritten as 
S 



dr 



8r 



U 

2-cos(6 J+1 -9 J ) 



•(7) 



If one takes the continuum limit of d — > and L — > 
oo while fixing the value of Ld, the action of Eq. ([7]) 
coincides with that for the spinless Tomonaga-Luttinger 
(TL) liquid gHESt 



S[9) 



hK 

27 



Ld ,M 



dx 



dr 



00 



86 

o.v 



.(8) 



where it is obvious that K is the TL parameter. Notice, 
however, that the values of K and u in Eq. ([5]) are renor- 
malized due to the effects of high-energy modes and Umk- 
lapp scattering and that the original relations of those pa- 
rameters with U / J and v no longer hold. The TL liquid 
model generally describes low-energy physics of a mass- 
less one-dimensional system. This indicates that while 
we analyze the discrete quantum rotor model of Eq. (J7JI , 
low-energy properties found in our analyses should be 
general in the spinless TL liquid. 

It is convenient to express the imaginary time in units 
of the Josephson plasma time h/Ej as 



T = T/> 



(9) 



where _Ej = hu/(y/2d) is the Josephson plasma energy. 
Inserting Eq. (j9]) into Eq. (J7J, we obtain 



S = h _ s, 
y/2ir 

where s is the dimensionless action 



df 



1 86 86 
28f ' df 



V(6) 



(10) 



(11) 



In order to express the action compactly, we introduced 
in Eq. ([TT]) an L-dimensional vector 9 defined as 



h(f ),..., 6 (f ),..., 6 L (f)f 



(12) 



and the potential, 



V0) = Y.VM+uOj) 

3 = 1 
L 

= ^-2cos(0 J+1 -^ 



(13) 



and /3 = /3Ej. From Eqs. (TT]) and (|10p we clearly see that 
h c = \/2tt I K plays the role of the effective dimensionless 
Planck's constant for this problem. The limit of h e — > 
corresponds to the classical (Gross-Pitaevskii) regime, 
while at h c > 1 quantum fluctuations become significant 
and can even drive the system to a different insulating 
phase. 



III. QUANTUM NUCLEATION RATE FOR THE 
PHASE SLIPS 



Extremizing the action by imposing 5s = 0, we obtain 
the classical equations of motion for the phases 6j , 



-2sin(0j+i - 6j) +2sin(0j - 0j_i) . (14) 



8 2 6 ] 
df 2 



There are two types of stationary solution of Eq. (JTj 
The first one describing a state carrying a homogeneous 
superflow with the winding-number n is 



(15) 



This state possesses the (quasi-)momentum p = 
2irhn/ (Ld) . The other is a saddle-point solution with a 
phase kink separating two states with different winding 
numbers: 



_ a, 
'n,j - 2 



l +^n(i-l), 



where 



and 



L - 1 + 2n 

a r , = — 7r mod 27r, 

L — 2 



2un — a r 
L - 1 



(16) 



(17) 



(18) 



The value of <p n defines the momentum in the system and 
a n is the phase difference between the 1st and L-th sites, 
the location of the phase kink. The magnitude of this 
kink a n is defined within the interval [— 2ir,0]. In the 
limit of the large number of sites L 3> n the expression 
for a n simplifies: 



-tt 1 + 



1 + 2n 
L 



mod 2tt. 



(19) 



For small winding numbers n <C L the phase kink is 
approximately equal to tt. 

We consider that a metastable flowing state with the 
winding number n is prepared at the initial time t = 0. 
The flow momentum is assumed to be smaller than the 
critical value p c = Hit/ (2d), above which the uniform 
flowing solutions become unstable. In the classical limit 
(h c — > 0) and at zero temperature the system remains 
in the initial state for an infinitely long time, i.e. the 
superflow is persistent. In contrast, when quantum fluc- 
tuations are strong enough, the metastable state decays 
into states with smaller momenta through quantum nu- 
cleation of phase slips and the lifetime of the superflow 
is finite. The decay rate of the metastable state, i.e. the 
nucleation rate of the quantum phase slip, can be calcu- 
lated using the celebrated instanton formula [43l - |47| : 



hT ~ E } LA 



SB 



SB 

exp [ — — 
2iTho \ hr 



(20) 
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Here sb is the action for the bounce solution 6(f) 
9 B (f ) of Eq. (fT4"]) , and the prefactor A is given by 



A = 



EL a 



(0) 



1/2 



(21) 



where A TO 's and Am 's are the solutions of the following 
eigenvalue equations: 



M£ m (f) = A m £ m (f), 



and 



(22) 
(23) 



^ (0) 4°Ht) = A»(t). 
Here the L-dimensional vectors 

U = (a,m(r), . • .,^-,m(f), . . , (24) 

and 

^ > = (f$S.(r),-..,CjS.(r),...,fgL(r)) t . (25) 
obey the orthonormalization condition 

dr 6 • f ro = fy m , / df $ 0) ■ = 6 l>m . (26) 



The L x L matrices M. and A^ ) are defined as 




d 2 d 2 V, 



df 2 d6 2 



d6 2 



(27) 



and 



Note that the factor of L in the 



1=82? 



where to 2 = da.V-i 



right hand side of Eq. (|20|) reflects the fact that there are 
L independent trajectories corresponding to the phase 
slip happening at one out of L links [3^, |4J] . We empha- 
size that the use of the quantum rotor model is advan- 
tageous in the sense that Sb and A do not depend on u 
and U / J but depend only on p and L. This means that 
frr/Ej depends on U, J, and v only through h c . 

To obtain the bounce action §b and the prefactor A, 
we numerically calculate the bounce solution B (f). The 
bounce solution 9 B (f) starts with the metastable state 
with the winding number n, i.e. 9 B (—00) = 6„ , goes 
through the saddle point 0%, bounces at the classical 
turning point, and returns to the initial state. An ex- 
ample of the bounce solution is shown in Fig. [T] for n = 1 
and L = 60, where the origin of the time is set such that 
the bounce solution reaches the classical turning point 
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FIG. 1: (Color online) (a) The bounce solution for the initial 
state with the winding number n = 1. We take L = 60. It 
is clear that the bounce solution forms a pair of a vortex and 
an anti- vortex in the (x, r)-plane. (b) The bounce solution is 
shown in section along the lines f — —30 (black circles), —3.75 
(red squares), and (blue diamonds), which correspond to 
the metastable state, the vicinity of the saddle point, and the 
classical turning point, (c) The bounce solution is shown in 
section along the line j = L/2. 



at f = 0. In Fig. Ufa), we see that the bounce solu- 
tion for the quantum phase slips forms a pair of a vortex 
and an anti- vortex in the (x, r)-plane. As clearly seen 
in Fig. Hlb), the phase kink is not located at a lattice 
site, but at a link between two sites (30th and 31st), and 
thereby the density rij at any sites does not vanish in 
contrast to the phase-slips in continuous systems. This 
is consistent with the basic assumption of the quantum 
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FIG. 2: (Color online) (a) The bounce action Sb versus the 
system size L. (b) The prefactor A versus the system size 
L. The (quasi-)momentum is fixed to be pd/h = 2n/5. The 
red solid line represents the best fit to the data with a fitting 
function f(x) — a + bx~ c . The blue dashed line represents 
the extrapolated value a. 

rotor model that the density fluctuation is small com- 
pared to the density itself. 

Substituting the bounce solution into Eq. we ob- 
tain the bounce action. Moreover, solving the eigenvalue 
equations Eqs. ([22]) and (|23|) with the bounce solution, 

we obtain £ OT and fJvi • Once £ m and §n are obtained, 
we can calculate the prefactor A through Eq. (f2~Tj). 

When calculating the bounce solution, one often intro- 
duces collective variables to reduce the number of degrees 
of freedom [U [22j. However, the use of such collec- 
tive variables restricts the analyses to a small region with 
respect to the momentum. In contrast, we deal with all 
the phase degrees of freedom, and this unbiased treat- 
ment enables us to more accurately obtain the nucleation 
rate for the entire region of the momentum. 

Let us calculate the bounce action sb and the prefactor 
A as functions of the momentum p. For a given value 
of p, sb and A depend on the number of lattice sites 
L. A typical example is shown in Fig. [2j where sb and 
A at pd/h = 2n/5 are plotted by varying L. When L 
increases, these two quantities are converged to certain 
asymptotic values. We extract the asymptotic values by 
fitting the numerical data to a function f(x) = a+bx~ c as 
represented by the blue lines in Fig.[3J This way allows us 
to obtain the values of sb and A for the thermodynamic 
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FIG. 3: (Color online) (a) The bounce action Sb as a function 
of the (quasi-)momentum p. The red squares represent Sb 
obtained by the extrapolation shown in Fig. [2] while the blue 
circles represent Sb for the system size L — 2nh/(pd), where 
the winding number is n — 1. The red solid line represents 
sb = 7.26(7r/2 — pd/K) 5 ^ 2 . The blue dashed line represents 
the best fit to the data for pd/h < 7r/18 with a function 
f(x) — alog(a:) + b, where the fitting parameters turn out to 
be a = -9.04 ± 0.05 and b = 3.84 ± 0.04. (b) Magnified view 
of (a) that focuses on the region near the critical momentum 
pd/h= 7r/2. 



limit (L — > oo). In Figs. [3] and HI the bounce action 
sb and the prefactor A versus the momentum p for the 
thermodynamic limit are plotted by the red squares. 

While the extrapolation to the thermodynamic limit 
is applicable for a region of relatively large momenta 
(pd/h > 7r/10 in Figs. [3] and H]), it is practically difficult 
for very small momenta because the calculations for very 
large systems are required. For this reason, in the region 
of small momenta we calculate sb and A only by taking 
the system size L = 2Trh/(pd), i.e. the winding number 
n = 1. For instance, we take L = 40 for pd/h = ir/20. 
In Figs. [3] and [4] §b and A for n = 1 are plotted by the 
blue circles. Although the values of Sb and A are a little 
overestimated without the extrapolation, the system size 
for a small momentum is so large that the deviation from 
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FIG. 4: (Color online) The prefactor A as a function of the 
(quasi-)momentum p is plotted in a log-log scale. The red 
squares represent A obtained by the extrapolation shown in 
Fig. [21 while the blue circles represent A for the system size 
L = 2nh/(pd), where the winding number is n = 1. The 
red solid line represents A = 0.867(7r/2 - pd/h) l/2 . The blue 
dashed line represents the best fit to the data for pd/h < 
7r/18 with a fitting function fix) — ax b , where the fitting 
parameters turn out to be (a,b) = (1.68±0.01, -1.87±0.01). 



the values in the thermodynamic limit is fairly small as 
already seen in the data points at pd/h = n/10 (L = 20) 
in Figs. [3] and |U 

In the vicinity of the critical momentum p c = Hit/ (2d), 
it has been predicted by previous work [22[ that the 
bounce action scales as Sb = C s (tt/2 — pd/h) 5 / 2 . In a 
similar way, the scaling formula for the prefactor can be 
obtained as A = CU(7r/2— pd/fi,) 1 / 2 . Inserting the numer- 
ical data of Sb and A for pd/H = 8ir/ 17 into these formu- 
lae, we obtain the coefficients C s — 7.26 and Ca = 0.867. 
As shown in Figs. [3J and HI the formulae with these val- 
ues of the coefficients agree well with the numerical data 
(red squares) for the momentum close to p c . Surprisingly, 
we find that the agreement in sb is almost perfect until 
pd/h ~ 7r/6 which is far away from p c . Note that while 
the previous work has predicted C s = 7.1 that is indeed 
close to our prediction, it is a little less accurate than 
ours because of the use of a variational ansatz [22j . 

For small momenta, p h/d, we numerically find that 
the bounce action exhibits a logarithmic dependence as 

s B =a s log(pd/h) + b s , (29) 

and that the prefactor obeys a power law as 

A = a A (pd/h) bA . (30) 

We extract the coefficients a s = —9.04, b s = 3.84, a A = 
1.68, and b A = —1.87 by fitting these formulae to the 
numerical data of sb and A for n = 1 represented by 
blue circles in Figs. [3] and HI Since the formula is valid 
for small momenta, we used the data in the region of 
pd/h < n/18 for the fitting. Substituting Eqs (|2T))) and 



(|30|) into Eq. ((20]) , the scaling formula for the decay rate 
is derived as hT/(LEj) cx (pd/h) 2mK ~ 1S7 . From this 
numerical result, we argue that the nucleation rate obeys 
the following power law, 



hv 

LE~, 



= Cri P 4 



2K-2 



(31) 



where the coefficient Cr is independent of the momen- 
tum p. The same scaling formula has been derived for ID 
homogeneous superconductor with energy dissipation at 
the phase-slip centers [l5[ . We emphasize that this agree- 
ment is not trivial because the model used in Ref. [l5| is 
qualitatively different from ours in the sense that our 
model explicitly includes the lattice and does not in- 
clude energy dissipation. As explained in Sec. HV] this 
agreement is rooted in the fact that both models exhibit 
the Berezinskii-Kosteritz-Thouless (BKT) transition at 
K = 2. 

Since the resistance R is related to the nucleation rate 
as R oc Y/p [l3j ]. Eq. (piTj) indicates that a ID Bose 
fluid can flow with almost no resistance in a periodic 
potential as long as the flow velocity is sufficiently small. 
This is consistent with the previous numerical result that 
the dipole oscillations of ID lattice bosons confined in a 
parabolic potential are hardly damped when the flow ve- 
locity is more than one order of magnitude smaller than 
the mean-field critical velocity [24[ ■ 

Note that in deriving Eq. (|3"Tj) we ignored the loga- 
rithmic correction stemming from the factor ^/sb. The 
formula of Eq. (|3"Tj) has been found through the numerical 
analyses on the discrete quantum rotor model Eq. (|10l) . 
However, it is very likely that this formula is generally 
valid for the spinless TL liquid in a periodic potential 
because it is a low-energy property. 



IV. QUALITATIVE DISCUSSIONS ON THE 
SCALING FORMULA 



In this section, we qualitatively explain a reason 
why the decay rate should obey the scaling formula of 
Eq. (fJTj) . Our explanation is twofold. First, the term 
2K in the exponent 2K — 2 can be understood as the 
contribution from the bounce action, whose analytical 
expression can be obtained by using the analogy with 
the classical 2D XY model. Secondly, the term —2 in the 
exponent can be determined by considering the relation 
between the quantum nucleation of phase slips and the 
quantum phase transition between the superfluid and the 
Mott insulator. 

Let us explain the first item. It is well known that 
taking the continuum limit and regarding ut as an- 
other spatial variable y, the ID quantum rotor action 
of Eq. (fTO)) is equivalent to the energy of the classi- 
cal 2D XY model. In this mapping, the bounce solu- 
tion corresponds to a vortex-antivortex pair. This anal- 
ogy allows one to express the bounce action in terms 
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FIG. 5: The size of a vortex-antivortex pair tb versus the 
inverse (quasi-)momentum H/(pd). The blue line represents 
the best fit to the data for pd/h < n/18 with a function 
f(x) — ax + b, where the fitting parameters turn out to be 
a — 0.753 and b = 0.412. Notice that the definition of tb is 
depicted in Fig. QJb). 



of the distance between the vortex and antivortex tb- 
When tb is sufficiently large compared to the size of 
vortex cores, the bounce action is well approximated as 
S^e/h = 2K log tb + s v , where s v is the contribution from 
the vortex cores and independent of tb [48]. In order 
to find the p-dcpendcncc of Sb , we numerically calculate 
tb versus the inverse momentum K/(pd) for the winding 
number n = 1 as shown in Fig. [5j There we observe that 
when pd/h <^ 1, tb linearly increases with h/(pd). Using 
this relation, we obtain the p-dependence of the bounce 
action as Sq/H = —2K\og(pd/h) + Const, and thereby 



(pd/h) 



2 A 



Hence, it is natural to anticipate 



r oc L(pd/h) 2K+x , where X is a constant that will be 
determined to be X = —2 in the following. 

In order to corroborate X = —2, we focus on the nu- 
clcation rate from the state with a certain fixed wind- 
ing number n, which is given by r„_j.„_i oc £ I - 2K — X + 1 . 
We discuss the relation between the nucleation rate and 
the quantum superfluid-insulator transition. It is impor- 
tant to remind us that the instanton formula is derived 
within the dilute gas approximation (DGA), in which 
the bounces are assumed to be well separated from each 
other [43M47| . This means that when DGA breaks down, 
many vortices are created in the space-time coordinate so 
that they destroy the long-range phase coherence, leading 
to the quantum phase transition to an insulating phase. 
In other words, the breakdown of DGA signals the Mott 
transition J4i| . In general, DGA is valid when the size of 
a bounce in the imaginary-time axis tb is much smaller 
than the nucleation time 1/T 0, In the present 

case, the bounce time is inversely proportional to the 
momentum as shown in Fig.[3J which means tb oc L, and 
thereby TBr„_j.„_i oc L~- 2K - X + 2 . This means that when 
—2K — X + 2 < 0, DGA is valid in the thermodynamic 



limit, i.e. the system is in the superfluid phase. There- 
fore, the condition — 2K — X + 2 = has to be fulfilled at 
the Mott transition point. From a different point of view, 
the renormalization group analyses have shown that the 
transition of the BKT type occurs at K = 2 in the pres- 
ence of a periodic potential [U |42[. Thus, we reach the 
conclusion that X = — 2. Notice that X = — 2 can be 
derived in the same way also in the model of diffusive 
ID superconductors studied in Ref. [TH from the fact that 
there exists the BKT transition at K = 2 as well. 

The explanation based on the relation between DGA 
and the superfluid-insulator transition is applicable also 
to the cases of a disorder potential and a single bar- 
rier. For the disorder case, Khlcbnikov and Pryadko have 
found that the nucleation rate scales asT oc L(p/p c ) 2K ^ 1 . 
Anticipating T oc L(p/p c ) 2K+x , let us corroborate that 
X = — 1 along the procedure introduced above. Since 
tb oc L as in the case of a periodic potential, the onset of 
the DGA breakdown, i.e. the transition point to an insu- 
lating phase, is given by —2K — X + 2 = 0. Meanwhile, 
it is known from renormalization group analyses that 
K = 3/2 at the transition to the Bose glass phase [50| . 
Therefore, a simple algebra leads to X — — 1. For the 
single- barrier case, it has been shown in Refs. [ill. fl9l that 
r oc {p/p c ) 2K ~ 1 ■ Notice the absence of the factor of L, 
which reflects the fact that the phase slips occur only at 
the barrier potential. Anticipating T oc [p/p c ) 2K+x , one 
can easily show that X = — 1 as follows. The onset of 
the DGA breakdown is given by — 2K — X + 1 = while 
renormalization group analyses have shown that a tran- 
sition to an insulating phase pinned by the barrier occurs 
at K = 1 [IJ. Hence, X = — 1. 

Furthermore, we briefly discuss the effects of finite 
temperatures on our results. For the scaling formula 
of Eq. (j5Tj) to be valid, the bounce time tb has to be 
sufficiently small compared to the inverse temperature 
h/{k-B,T), Since tb ~ h/(pd) as shown in Fig. [SI this con- 
dition turns out to be k^T Ejpd/h. In the interme- 
diate temperature region, namely Ejpd/h <C k^,T <C Ej, 
the bounce time is bounded by the inverse temperature 
as tb ~ Ej/(k-aT). In this case, the nucleation of phase 
slips occurs due to the thermally assisted quantum tun- 
neling [52l | , and the nucleation rate is given by 



r oc l- 



pd ( ksT 

~e7 



2K- 



(32) 



where the constant c has been determined to be c = 
3 in previous work [53[. Notice that scaling formulae 
similar to Eq. (|52l have been found also in the cases of a 



disorder potential [21[ and a single-barrier potential |18l . 
EU. Equation ([32]) means that the transport is ohmic, 
but that the resistance can be very small when K 3> 1. 
This ensures the presence of superfluidity in the practical 
sense even at small finite temperatures When the 

temperature is as high as k^T S> £j, the phase slips 
can not be nucleated in the space-time coordinate and 
the thermal activation process becomes dominant to the 
supcrflow decay. 
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FIG. 6: (Color online) The red solid line represents the time 
evolution of the flow velocity v(t) in the dynamics of the ID 
BHM, where L = 160, v = 1, U/J = 3, and n = 4. The 
blue dashed line represents the flow velocity at the winding 
number n = 3. 



,1.00 




FIG. 7: (Color online) The red solid line represents the time 
evolution of the persistence probability P(t) in the dynamics 
of the ID BHM, where L = 160, v = 1, U/J = 3, and n = 4. 
The longitudinal axis is shown in a logarithmic scale. In the 
region sandwiched between the two green dotted lines, P(t) 
decays exponentially. The blue dashed line represents the best 
fit with a function of Eq. 1)35[) to the data in the exponentially 
decaying region. 



Since the Josephson plasma energy Ej separates the 
quantum-tunneling regime from the thermal-activation 
regime, it is important to provide an estimate of Ej for 
present typical experiments. Taking the experiment of 
Ref. for example, Ej/ks ~ 30nK because the sound 
velocity is u ~ 2.1mm/s and the lattice spacing is 
d = 405nm. Given the fact that recent experiments have 
lowered the temperature as low as 5nK [3(1 [54[ , it is very 
likely that the regime of /cbT <C Ej is experimentally ac- 
cessible so that one can observe the superflow decay via 
quantum nuclcation of phase slips. 



V. TEBD ANALYSES OF THE 
BOSE-HUBBARD MODEL 



potential with an integer filling for the following two rea- 
sons. First, this formula is a property at small momenta 
(jp -C p c ), and such a low-energy property should be valid 
commonly in the spinless TL liquid. For example, the 
compressibility and the long-range behaviors of correla- 
tion functions are expressed in terms of the TL parame- 
ter K and the sound velocity u regardless of microscopic 
details of the original Hamiltonian [4Lj] . Secondly, as dis- 
cussed in the previous section, this formula is closely re- 
lated to the Mott transition at K = 2, which is a common 
property in the spinless TL liquids with an integer filling 
in a periodic potential. 

In this section, we study the superflow decay in the 
Bosc-Hubbard model with unit filling in order to demon- 
strate that the applicability of the scaling formula is 
not limited to the quantum rotor model. Recall that 
the quantum rotor model quantitatively agrees with the 
Bose-Hubbard model only in the region of high filling 
factors (y ^> 1) [Hj]. It is important to investigate the 
unit-filling case also because the supcrfluid transport of 
ID Bose gases in optical lattices has been experimen- 
tally studied mainly in the region of low-filling factors 
(y ~ 1). In the low-filling regime, it has been predicted 
within the GP mean field theory that the Landau insta- 
bility, which is characterized by the emergence of excita- 
tions with negative energies, sets in at momenta smaller 
than the critical value for the dynamical instability, i.e. 
p = Kit/ (2d) [5^, [5(|. However, this instability can break 
down superflow only when the temperature is sufficiently 
high so that the thermal fraction is comparable to the 
condensate fraction (57l - [59| . and is not relevant in the 
regime of our interest in which quantum fluctuations are 
dominant over thermal ones. 

We describe ID lattice bosons in a ring-shaped geom- 
etry with the following BHM with a phase twist: 



H= -J 



— ^hjihj - 1),(33) 

" 3=1 



While the scaling formula of Eq. (|3Tj) was obtained 
from the 0(2) quantum rotor model, it is expected to 
hold generally for ID spinless supcrfluids in a periodic 



where bL+i = b±, reflecting the nature of the ring-shaped 
geometry. In Eq. (|33[) . we include the phase twist e~ %B in 
the hopping term in order to control the winding num- 
ber of states. To deal with the quantum dynamics of 
superflow of the ID BHM Eq. ([33]), we use the TEBD 
method [33| for a periodic boundary condition [34| , which 
allows us to accurately compute the time evolution of 
many-body wave functions in ID quantum lattice sys- 
tems. It has been shown in our previous work that TEBD 
is applicable to the problem of superflow dynamics associ- 
ated with quantum phase slips [39|] . We first calculate the 
ground state of Eq. (|33| with the phase twist 9 = 2-nn/L 
via the imaginary time evolution, and thereby a flowing 
state with the winding number n that is metastable in 
the classical limit is prepared. Taking this state as the 
initial state and setting 9 = at t = 0, we compute the 
real-time evolution. 

In Fig. |51 wc show the time evolution of the averaged 
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FIG. 8: The red circles represent the nucleation rates of quan- 
tum phase slips V extracted from the real-time dynamics of 
the ID Bose-Hubbard model with L — 160 as functions of the 
flow (quasi-)momentum p, where U/J — 2.8 (a), 3 (b), and 
3.2 (c). The plots are shown in a log-log scale. In each plot, 
the blue solid line represents the scaling formula of Eq. (|3ip 
with the constant Cr determined such that the line passes on 
the data point with the smallest momentum. The TL param- 
eters are taken from Ref. l6ll as K = 2.52 (a), 2.37 (b), and 
2.17 (c). 



flow velocity is constant in time if one computes classi- 
cal dynamics of the Gross-Pitaevskii equation neglecting 
quantum fluctuations. 

To quantify the tunneling rate from the metastable 
state, i.e. the nucleation rate of quantum phase slips, we 
calculate the overlap of the wave function with the initial 
state P{t) = |($(<)|^(0))| 2 , which can be interpreted as 
the persistence probability, i.e., the probability that the 
state remains in the initial state by the time t. It is ex- 
pected from the tunneling theory that when the initial 
state is a metastable state, the persistence probability 
decays exponentially as P(t) ~ exp(— Yt) with the nu- 
cleation rate T [(3(3. Notice that the ability to calculate 
the persistence probability is a clear advantage of TEBD 
for a periodic boundary condition over TEBD for infinite 
systems that was used in Ref. I23I . 

In Fig. [3 we show P(t) with the same parameters as 
used in Fig. [6] Indeed there is a large region where P(t) 
exhibits the exponential decay Wc extract the nucleation 
rate Y by fitting the data in the exponentially decaying 
region to the following function, 



g(t) = Dcxp(-Yt), 



(35) 



and taking D and Y as free parameters. In Fig. [HJ wc plot 
the nucleation rates versus the momentum p for three 
values of U / J, where L = 160. We also show the scaling 
formula of Eq. (f3"Tj) represented by the (blue) solid lines, 
where Cr is taken such that the line passes on the data 
point with the smallest momentum. We use the TL pa- 
rameter K numerically extracted from the single-particle 
correlation function (6j6o) with the distance 16 < r < 32 
in Ref. l6l|. Notice that the TL parameter K in the 
present paper is equivalent to 1/K in the definition of 
Ref. [HH. Wc clearly sec that the data points approach 
the lines when the momentum decreases, thus justifying 
the validity of the scaling formula for BHM with unit 
filling. 



VI. CONCLUSIONS 



flow velocity, 

" = |*-^ (34) 

3 

for L = 160, U/J = 3, and n = 4. We see that the 
flow velocity decreases in time, clearly exhibiting the su- 
perflow decay due to quantum tunneling. However, the 
averaged flow velocity docs not exhibit a sudden drop 
by a quantized amount, which could be a characteristic 
of phase slips, but gradually decreases in time. This is 
because the phase slip jump is smoothened out by tak- 
ing the quantum average of many events. In each event 
a phase slip occurs at a different time. Notice that the 



In summary, we have studied the decay of supcrflow 
via quantum nucleation of phase slips in onc-dimcnsional 
(ID) superfluids in the presence of a periodic potential. 
Within the quantum rotor regime, we used the instanton 
method to obtain the nucleation rate for all the region 
of the momentum p. When the momentum is close to 
the mean-field critical value p Cl we improved the expres- 
sion of nucleation rate that was previously obtained in 
Ref. I22I . For small momenta p <C p c , we derived the 
scaling formula of the nucleation rate with respect to p, 
which is expressed in Eq. (|3"Tj) . We discussed the relation 
between the dilute gas approximation and the quantum 
superfluid-insulator transition in order to gain a unified 
physical interpretation of the scaling formulae for peri- 
odic, disorder, and single-barrier potentials. Applying 
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the time-evolving block decimation method to the ID 
Bose-Hubbard model with unit filling, we analyzed the 
quantum dynamics of the supcrflow decay and confirmed 
the validity of the scaling formula of Eq. (j3"Tj) . 

While we have calculated the nucleation rate of quan- 
tum phase slips in order to characterize the supcrflow 
decay, it still remains ambiguous how the nucleation rate 
is related to the transport of ID Bose gases in the pres- 
ence of a trapping potential that has been studied in cold 
atom experiments [8l4ll| . Since the damping rate of the 
dipole oscillations has been often used to quantify the 
transport of trapped atomic gases, in our future work we 
will clarify direct connections of the nucleation rate with 



the damping rate. 
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